bamstats

first we look at coverage of the data going into stacks, the bam files

pulled up the bamstats from log output of gstacks (“gstacks.log.distribs”), then edited file with regex include a population flag

bamstats <- read_tsv("log_files/12-8_run/12-8bamstats.txt")
ggplot(data = bamstats)+geom_density(aes(x=primary_kept))+geom_vline(aes(xintercept = median(bamstats$primary_kept)), color = "red")+ggtitle("Retained record from BAM input")

bamstats %>%
  group_by(pop) %>%
  summarise(median_forward_reads = median(primary_kept), mean_kept = mean(kept_frac) )

Here we see no issue with variation in mapping quality across the populations and good evidence that we got more than the number of reads we were looking for per individual. ~13 million total reads per samples (or 7.5million pairs), which if we have about 350k rad loci works out to about 20x coverage before filtering sites (but after filtering for bad reads in the process_radtags fork).

We lost about 20% of the reads in the bam files to default quality filters in the gstacks fork (minmapq 10, max-clipped 0.20, unmapped reads) summary below:

Read 6556590170 BAM records:
kept 5322753707 primary alignments (83.4%), of which 2628057004 reverse reads
skipped 423599169 primary alignments with insufficient mapping qualities (6.6%)
skipped 455204465 excessively soft-clipped primary alignments (7.1%)
skipped 179611285 unmapped reads (2.8%)
skipped some suboptimal (secondary/supplementary) alignment records

At this point the expected coverage is ~ 17x (2433922417 paired reads per 395 samples (about 6 million reads per sample) per estimated 350k rad loci)

Raw Genotypes Coverage

similar to above we pulled the distribution of effective coverage from the stacks logs, edited to include a population flag and examined below

First, the summary output from gstacks.log:

Built 1000251 loci comprising 2694696703 forward reads and 2433922417 matching paired-end reads; mean insert length was 244.1 (sd: 75.2).
Removed 260774286 unpaired (forward) reads (9.7%); kept 2433922417 read pairs in 962035 loci.
Removed 796945198 read pairs whose insert length had already been seen in the same sample as putative PCR duplicates (32.7%); kept 1636977219 read pairs.

Genotyped 962035 loci:
effective per-sample coverage: mean=9.4x, stdev=5.1x, min=1.0x, max=36.9x
mean number of sites per locus: 523.7
a consistent phasing was found for 36695500 of out 38359065 (95.7%) diploid loci needing phasing

Quick summary: 33% PCR duplicates + vastly more rad loci than expected leads to drastic reduction in coverage… but lets look more closely

eff_cov <- read_tsv("log_files/12-8_run/effective_cov_12-8.txt")
a <- ggplot(data=eff_cov)+geom_density(aes(x= n_loci))+geom_vline(aes(xintercept=median(eff_cov$n_loci)), color = "red")+ggtitle("Number of Rad Loci")
b <- ggplot(data=eff_cov)+geom_density(aes(x= mean_cov_ns))+geom_vline(aes(xintercept=median(eff_cov$mean_cov_ns)), color = "red")+ggtitle("weighted mean coverage")
c <- ggplot(data=eff_cov)+geom_density(aes(x= pcr_dupl_rate))+geom_vline(aes(xintercept=median(eff_cov$pcr_dupl_rate)), color = "red")+ggtitle("pcr duplicates")
d <- ggplot(data = eff_cov)+geom_point(aes(pcr_dupl_rate, mean_cov_ns), alpha = 0.1) +geom_smooth(aes(pcr_dupl_rate, mean_cov_ns), method = "lm") 
plot_grid(a,b,c,d)

rm(a,b,c,d)

At this point, things are not looking great, about 50% of samples have less than 10x effective coverage, and this seems largely to be caused by (1) a great increase in the number of rad loci realtive to what was expected given other papers numbers for the number of PstI cut sites and my own in silico digest and (2) PCR duplicates take up anywhere from ~18% to 60% of the reads depending on the library (it’s clear that the variation in pcr duplication rate is explained by variation at the level of library prep given the clustering in the last plot above).

quick note: a previous examination of PCR duplication reveals that only 0.6% of duplicates appear to be optical duplicates, suggesting that library complexity is the root issue here, not sequencing.

The good news:
if there are actually >500k PstI cut sites in the stickleback genome, than we have more sites than we need to effectively cover the whole genome - In the worst case from the literature (Roesti et al 2015), LD decays at 1kb, and the stickleback genome is 447mb, so achieving full genomewide coverage requires 450k radloci (IN THE WORST POSSIBLE CASE), therefore we can probably lose some of the poorly sequenced sites and individuals to improve depth and still have markers that cover the entire genome!

To look at things at this level we will have to examine coverage per site per individual, instead of these summary statistics, this is below

r0.95

this section covers the first fun of populations (-r 0.95 i.e. only retain variable sites if rad locus is covered in in 95% of sample within a population)

How many of the rad loci are “bad” (i.e. absent int >5% of samples within a population)?
>Removed 653042 loci that did not pass sample/population constraints from 962035 loci.
Kept 308993 loci, composed of 248618752 sites; 2268754 of those sites were filtered, 838333 variant sites remained.
221216030 genomic sites, of which 26716049 were covered by multiple loci (12.1%). Mean genotyped sites per locus: 795.32bp (stderr 0.28).

How does coverage look now?

#first capture the depth per locus per individual using vcftools

#!/bin/bash

# set max wall-clock time (D-HH:MM:SS)
#SBATCH --time=0-23:59:00

#SBATCH --cpus-per-task=1

# set partition/queue to use
#SBATCH --partition=day-long-cpu

# set name of job
#SBATCH --job-name=depths

# set name of output file
#SBATCH --output=depths.out

/opt/vcftools_0.1.13/vcftools  --vcf /home/ddayan/plastic/stacks/genotypes/populations.snps.vcf --depth 

/opt/vcftools_0.1.13/vcftools  --vcf /home/ddayan/plastic/stacks/genotypes/populations.snps.vcf --site-depth 

/opt/vcftools_0.1.13/vcftools  --vcf /home/ddayan/plastic/stacks/genotypes/populations.snps.vcf --site-mean-depth 

/opt/vcftools_0.1.13/vcftools  --vcf /home/ddayan/plastic/stacks/genotypes/populations.snps.vcf --geno-depth 

#then randomly sample this geno-depth output to a reasonable size for plotting (50k snps)
shuf -n 10000 out.gdepth > gdepth_10k
( sed 1q out.gdepth; sed 1d gdepth_10k ) > gdepth_10k.txt 

gdepths <- read_tsv("log_files/12-8_run/gdepth_10k.txt")
idepth <- read_tsv("log_files/12-8_run/out.idepth")
ldepthmean <- read_tsv("log_files/12-8_run/out.ldepth.mean")
gs <- as.data.frame(unlist(gdepths[c(3:397)]))
colnames(gs) <- "depth"

a <- ggplot(data=idepth)+geom_density(aes(x=MEAN_DEPTH))+geom_vline(aes(xintercept = median(MEAN_DEPTH)))+ggtitle("Mean Depth per Locus per Individual")
b <- ggplot(data=ldepthmean)+geom_density(aes(x=MEAN_DEPTH))+geom_vline(aes(xintercept = median(MEAN_DEPTH)))+ggtitle("Mean Depth per Locus")+scale_x_continuous(breaks = scales::pretty_breaks(n = 10))
c <- ggplot()+geom_density(aes(x=gs$depth))+scale_x_continuous(breaks = scales::pretty_breaks(n = 10))+xlim(0,50)

plot_grid(a,b,c)

Something not good is going on here, when we look at the SNP level more than 50% of the sites have > ~12x coverage, suggesting some coverage based filtering could produce a sizeable snp dataset to conduct our analyses on, but when we look at the SNPxIndividual level we see that many SNPs have 0 depth. Let’s look more closely…

#lets make a clustered heatmap of coverage split across populations

smallgs <- sample_n(gdepths, 1000)
# smallgs$ID <- seq.int(nrow(smallgs))
# 
# 
# long_gs <- gather(gdepths, individual, depth, bb_f_10:rs_s_87)
# long_gs <- long_gs[complete.cases(long_gs),]
# 
# small_long_gs <- gather(smallgs, individual, depth, bb_f_10:rs_s_87)
# small_long_gs <- small_long_gs[complete.cases(small_long_gs),]
# samples <- sample(colnames(gdepths[,3:397]), 40)
# small_long_gs <- subset(small_long_gs, individual %in% samples)
# small_long_gs$pop <- substr(small_long_gs$individual, start=1, stop = 2)
# 
# #now lets reorder the rows according to clustering
# hc <- hclust(dist(as.matrix(small_long_gs[,c("depth", "ID", "individual")])))
# ord <- order.dendrogram(as.dendrogram(hc))
# small_long_gs$ID <- factor(x = small_long_gs$ID, levels = unique(small_long_gs$ID[ord]), ordered = TRUE)
# 
# 
# ggplot(data=small_long_gs)+geom_tile(aes(x=individual,y=ID, fill=depth, group = pop))+scale_x_discrete(expand=c(0,0))+scale_y_discrete(expand=c(0,0))+scale_fill_viridis_c()+facet_grid(. ~ pop, scales = "free")

#all that was very nice, but it turns out theres a great package that will do a better job of that
require(heatmaply)
samples <- sample(colnames(gdepths[,3:397]), 200)
smallergs <- smallgs[,samples]
smallergs <- t(smallergs)
smallergs <- smallergs[order(row.names(smallergs)),]
smallergs <- t(smallergs)
heatmaply(smallergs, Colv=FALSE, col_side_colors = substr(colnames(smallergs),start=1, stop = 2), titleX = FALSE, titleY=FALSE)

The plot above is a heatmap of depth with SNPs on the y axis and individuals (grouped by population (BB, CL, LB, RS)) on the x axis for a random subset of 200 individuals and 1000 snps
* Few SNPs have good coverage across all populations * BB is missing many rad loci, but coverage is good where a locus is retained
* RS is also looking rough
* Good overlap between CL and LB relative to RS and LB, which is odd because LB is recently diverged from RS

Questions: (1) does coverage overall vary between populations?
(2) Is zero coverage due to the r 0.95 filter, and if so, does it make more sense to filter dataset-wide rather than population-wise because this will lead to population level structure in missing data?
(3) alternatively, are these rad loci not present some populations (i.e. allele dropout), this seems like a LOT of cut site loss, and the fact that LB and RS are not more similar doesn’t fit with this idea

Rather than spend time going down this road, lets run populations again with no population wise filtering and see what we get.

R95

This run reatined only rad loci present in 95% of samples (-R 0.95)

ran the same vcftools depth utils and subsampled the big one

idepth.1 <- read_tsv("log_files/12-8_run/R95/out.idepth")
ldepthmean.1 <- read_tsv("log_files/12-8_run/R95/out.ldepth.mean")
gdepths.1 <- read_tsv("log_files/12-8_run/R95/gdepth_1k.txt")
gs.1 <- as.data.frame(unlist(gdepths.1[c(3:397)]))
colnames(gs.1) <- "depth"

a <- ggplot(data=idepth.1)+geom_density(aes(x=MEAN_DEPTH))+geom_vline(aes(xintercept = median(MEAN_DEPTH)))+ggtitle("Mean Depth per Locus per Individual")
b <- ggplot(data=ldepthmean.1)+geom_density(aes(x=MEAN_DEPTH))+geom_vline(aes(xintercept = median(MEAN_DEPTH)))+ggtitle("Mean Depth per Locus")+scale_x_continuous(breaks = scales::pretty_breaks(n = 10))
c <- ggplot()+geom_density(aes(x=gs.1$depth))+scale_x_continuous(breaks = scales::pretty_breaks(n = 10))+xlim(0,50)+geom_vline(aes(xintercept = median(gs.1$depth)))+ggtitle("Depth per Locus per Individual")

plot_grid(a,b,c)

This looks about the same as the r95 run in the first two plots, but the difference becomes apparent when the SNP_X_individuals with zero coverage are filtered away. This is accomplished (see below) by getting rid of 42% of the Rad loci from the first run.

Removed 782350 loci that did not pass sample/population constraints from 962035 loci.
Kept 179685 loci, composed of 153251877 sites; 2037395 of those sites were filtered, 676463 variant sites remained.
149138121 genomic sites, of which 4063490 were covered by multiple loci (2.7%).
Mean genotyped sites per locus: 844.50bp (stderr 0.31).

Lets see what this looks like in the heat map

#lets make a clustered heatmap of coverage split across populations

smallgs <- sample_n(gdepths.1, 500)
require(heatmaply)
samples <- sample(colnames(gdepths.1[,3:397]), 200)
smallergs <- smallgs[,samples]
smallergs <- t(smallergs)
smallergs <- smallergs[order(row.names(smallergs)),]
smallergs <- t(smallergs)
heatmaply(smallergs, Colv=FALSE, col_side_colors = substr(colnames(smallergs),start=1, stop = 2), titleX = FALSE, titleY=FALSE)

Much Better! Let’s just check one more time that the coverage and number of retained sites is consistent across populations

idepth.1$pop <- substr(idepth.1$INDV, 1, 2)
a <- ggplot(data=idepth.1)+geom_violin(aes(pop,MEAN_DEPTH))
b <- ggplot(data=idepth.1)+geom_violin(aes(pop,N_SITES))
plot_grid(a,b)

rm(a,b)

3

# make the new popmap of good inds
bamstats %>%
  filter(primary_kept > mean(primary_kept)-(2*sd(primary_kept))) %>%
  select(sample, pop) %>%
  write_tsv("./log_files/12-8_run/12-8_good_inds_popmap.txt", col_names = FALSE)

then ran populations with no filtering and this popmap which output then follwoing:

Removed 0 loci that did not pass sample/population constraints from 962035 loci. Kept 962035 loci, composed of 515209893 sites; 116 of those sites were filtered, 7350383 variant sites remained.
325605580 genomic sites, of which 141745630 were covered by multiple loci (43.5%). Mean genotyped sites per locus: 523.39bp (stderr 0.27).

7.3 million SNPs across 325 million bp

4 low confidence SNPs

#!/bin/bash

# set max wall-clock time (D-HH:MM:SS)
#SBATCH --time=0-23:59:00

#SBATCH --cpus-per-task=1

# set partition/queue to use
#SBATCH --partition=day-long-cpu

# set name of job
#SBATCH --job-name=depths

# set name of output file
#SBATCH --output=depths.out

/opt/vcftools_0.1.13/vcftools  --vcf /home/ddayan/plastic/stacks/genotypes/no_filt/populations.snps.vcf --minDP 5 --min-meanDP 10 --recode-INFO-all --recode --out 0.2

This step retained 13% of SNPs
981005 out of a possible 7350383 Sites

5 iterative missingness filter

#output missingess
/opt/vcftools_0.1.13/vcftools  --vcf /home/ddayan/plastic/stacks/genotypes/filtering/0.2.recode.vcf  --missing-indv

#get inds with >50% missingess
awk '$5 > 0.9 && NR>1 {print $1}' out.imiss > bad_inds1

no inds with >90% missing data

#!/bin/bash

# set max wall-clock time (D-HH:MM:SS)
#SBATCH --time=0-23:59:00

#SBATCH --cpus-per-task=1

# set partition/queue to use
#SBATCH --partition=day-long-cpu

# set name of job
#SBATCH --job-name=iterative_missingness

# set name of output file
#SBATCH --output=iterative_missingness.out

#filter out bad inds then filter bad genos
/opt/vcftools_0.1.13/vcftools  --vcf /home/ddayan/plastic/stacks/genotypes/filtering/0.2.recode.vcf  --remove bad_inds1 --max-missing 0.5 --recode-INFO-all --recode --out 0.3

After filtering, kept 980166 out of a possible 981005 Sites

#output missingess
/opt/vcftools_0.1.13/vcftools  --vcf /home/ddayan/plastic/stacks/genotypes/filtering/0.3.recode.vcf  --missing-indv

#get inds with >50% missingess
awk '$5 > 0.5 && NR>1 {print $1}' out.imiss > bad_inds2

removed 10 inds (4 bb, 5 lb, 1 rs)

#!/bin/bash

# set max wall-clock time (D-HH:MM:SS)
#SBATCH --time=0-23:59:00

#SBATCH --cpus-per-task=1

# set partition/queue to use
#SBATCH --partition=day-long-cpu

# set name of job
#SBATCH --job-name=iterative_missingness

# set name of output file
#SBATCH --output=iterative_missingness.out

#filter out bad inds then filter bad genos
/opt/vcftools_0.1.13/vcftools  --vcf /home/ddayan/plastic/stacks/genotypes/filtering/0.3.recode.vcf  --remove bad_inds2 --max-missing 0.8 --recode-INFO-all --recode --out 0.4

After filtering, kept 891407 out of a possible 980166 Sites

How many rad loci?

cut -f 3 0.4.recode.vcf | awk -F : '{print $1}' |  uniq | wc -l 

152973 unique catlog call numbers

coverage and het miscall rate check

before we do any filtering to explicitly try to find miscalled hets due to paralogues and other sources of false orthology at an assembled rad locus, lets see what our coverage looks like now and assess the heterozygote miscall rate

#make sample vcf
head -15 /home/ddayan/plastic/stacks/genotypes/filtering/0.4.recode.vcf > 0.4sampled.vcf
shuf -n 1000 /home/ddayan/plastic/stacks/genotypes/filtering/0.4.recode.vcf >> 0.4sampled.vcf

#!/bin/bash

# set max wall-clock time (D-HH:MM:SS)
#SBATCH --time=0-23:59:00

#SBATCH --cpus-per-task=1

# set partition/queue to use
#SBATCH --partition=day-long-cpu

# set name of job
#SBATCH --job-name=depths

# set name of output file
#SBATCH --output=depths.out


/opt/vcftools_0.1.13/vcftools  --vcf /home/ddayan/plastic/stacks/genotypes/final_12-8_filt/0.7.sampled.vcf --depth 

/opt/vcftools_0.1.13/vcftools  --vcf /home/ddayan/plastic/stacks/genotypes/final_12-8_filt/0.7.sampled.vcf --geno-depth 

gdepths <- read_tsv("log_files/12-8_run/gdepth0.4_1k.txt")
idepth <- read_tsv("log_files/12-8_run/0.4.idepth")
gs <- as.data.frame(unlist(gdepths[c(3:370)]))
colnames(gs) <- "depth"

a <- ggplot(data=idepth)+geom_density(aes(x=MEAN_DEPTH))+geom_vline(aes(xintercept = median(MEAN_DEPTH)))+ggtitle("Mean Depth per site per Individual")+xlim(0,50)
c <- ggplot()+geom_density(aes(x=gs$depth))+scale_x_continuous(breaks = scales::pretty_breaks(n = 10))+xlim(0,80)+ggtitle("Depth Per site per individual")

plot_grid(a,c, ncol = 1)

vcf <- read.vcfR(file = "log_files/12-8_run/0.4sampled.vcf")

## Scanning file to determine attributes.
## File attributes:
##   meta lines: 14
##   header_line: 15
##   variant count: 1000
##   column count: 377
## 
Meta line 14 read in.
## All meta lines processed.
## gt matrix initialized.
## Character matrix gt created.
##   Character matrix gt rows: 1000
##   Character matrix gt cols: 377
##   skip: 0
##   nrows: 1000
##   row_num: 0
## 
Processed variant 1000
Processed variant: 1000
## All variants processed

#now split this into populations to avoid wahlund effect
vcf_popsub <- function(vcf, pop, in_pop) {
  ids <- which(pop %in% in_pop)
  ids <- ids
  vcf <- vcf[, c(1,ids)]

  return(vcf)

}
popid <- unlist(substr(colnames(vcf@gt), 1,2))
vcf_cl <- vcf_popsub(vcf, popid, "cl")
vcf_bb <- vcf_popsub(vcf, popid, "bb")
vcf_rs <- vcf_popsub(vcf, popid, "rs")
# first get compute expected and observed genotype frequencies
gfreqs_cl <- exp_and_obs_geno_freqs(vcf_cl)
gfreqs_bb <- exp_and_obs_geno_freqs(vcf_bb)
gfreqs_rs <- exp_and_obs_geno_freqs(vcf_rs)

# then plot those.  Set max_plot_loci so that all 2000
# loci will be plotted
geno_freqs_scatter(gfreqs_cl, max_plot_loci = 1000)

geno_freqs_scatter(gfreqs_bb, max_plot_loci = 1000)

geno_freqs_scatter(gfreqs_rs, max_plot_loci = 1000)

#infer het miscall rate
overall_cl <- infer_m(vcf_cl, minBin = 1e15)
overall_cl$m_posteriors

overall_bb <- infer_m(vcf_bb, minBin = 1e15)
overall_bb$m_posteriors

overall_rs <- infer_m(vcf_rs, minBin = 1e15)
overall_rs$m_posteriors

whoa, a het zygote miscall rate of ~2%, not bad

6 Paralogs

First we remove sites with an unexpected allele balance at heterozygotes using picard and GATK

# to run picard
/opt/java/openjdk-1.8.0/bin/java -jar /opt/picard/build/libs/picard.jar 

#setting it up as an alias
picard="/opt/java/openjdk-1.8.0/bin/java -jar /opt/picard/build/libs/picard.jar"
$picard

#filterVcf - this doesnt work yet used gatk below
$picard FilterVcf  I=/home/ddayan/plastic/stacks/genotypes/filtering/0.4.recode.vcf O=0.5.vcf MIN_AB=0.1 R=/home/ddayan/plastic/stacks/genome/Gasterosteus_aculeatus.BROADS1.dna_sm.toplevel.fa

#but needs a sequence dictionary? 
/opt/java/openjdk-1.8.0/bin/java -jar /opt/picard-tools-1.119/CreateSequenceDictionary.jar R= /home/ddayan/plastic/stacks/genome/Gasterosteus_aculeatus.BROADS1.dna_sm.toplevel.fa O= /home/ddayan/plastic/stacks/genome/Gasterosteus_aculeatus.BROADS1.dna_sm.toplevel.dict

#and it needs to be indexed by samtools (samtools faidx)

# ... and it needs to sort the Vcf file itself first ...
$picard SortVcf \
      I=/home/ddayan/plastic/stacks/genotypes/filtering/0.4.recode.vcf \
      SEQUENCE_DICTIONARY=/home/ddayan/plastic/stacks/genome/Gasterosteus_aculeatus.BROADS1.dna_sm.toplevel.dict \
      O=0.4.sorted.vcf
      

#now anotate the genotypes with the allele balance tag
GATK="/opt/java/openjdk-1.8.0/bin/java -jar /opt/Gatk/GenomeAnalysisTK.jar "
$GATK    -R /home/ddayan/plastic/stacks/genome/Gasterosteus_aculeatus.BROADS1.dna_sm.toplevel.fa \
  -V /home/ddayan/plastic/stacks/genotypes/filtering/0.4.sorted.vcf \
  -T VariantAnnotator \
  -o /home/ddayan/plastic/stacks/genotypes/filtering/0.4.annotated.vcf \
  --annotation AlleleBalance \
  -U ALLOW_SEQ_DICT_INCOMPATIBILITY

#filter 
$picard FilterVcf  I=/home/ddayan/plastic/stacks/genotypes/filtering/0.4.annotated.vcf O=0.5.vcf MIN_AB=0.2 

Next we remove the bad sites (AB < 0.2) and calculate the depths at remaining sites

#now that we have filter flags done, we need to remove the bad rows and prepare to finish paralog filtering by calculating the mean site depths 

#!/bin/bash

# set max wall-clock time (D-HH:MM:SS)
#SBATCH --time=0-23:59:00

#SBATCH --cpus-per-task=1

# set partition/queue to use
#SBATCH --partition=day-long-cpu

# set name of job
#SBATCH --job-name=depths

# set name of output file
#SBATCH --output=depths.out


/opt/vcftools_0.1.13/vcftools  --vcf /home/ddayan/plastic/stacks/genotypes/filtering/0.5.vcf --remove-filtered-all --recode-INFO-all --recode --out 0.6

#############################

#!/bin/bash

# set max wall-clock time (D-HH:MM:SS)
#SBATCH --time=0-23:59:00

#SBATCH --cpus-per-task=1

# set partition/queue to use
#SBATCH --partition=day-long-cpu

# set name of job
#SBATCH --job-name=depths

# set name of output file
#SBATCH --output=depths.out


/opt/vcftools_0.1.13/vcftools  --vcf /home/ddayan/plastic/stacks/genotypes/filtering/0.6.recode.vcf --site-depth

After filtering, kept 872741 out of a possible 891407 Sites (AB filter)

Here we filter out sites with 2.5x the modal depth (2095 sites)

ldepth0.6 <- read_tsv("log_files/12-8_run/0.6.ldepth")

getMode <- function(x) {
keys <- unique(x)
keys[which.max(tabulate(match(x, keys)))]
}


ggplot()+geom_density(data=ldepth0.6, aes(SUM_DEPTH)) +geom_vline(aes(xintercept=(2.5*getMode(ldepth0.6$SUM_DEPTH))), color = "red") +geom_vline(aes(xintercept=(getMode(ldepth0.6$SUM_DEPTH))), color = "blue") 

tmp <- ldepth0.6[ldepth0.6$SUM_DEPTH > 2.5*getMode(ldepth0.6$SUM_DEPTH),c(1,2)]
bad_sites <- tmp %>%
  unite(col = "id", CHROM, POS, sep = "_")
write_tsv(bad_sites, "./log_files/12-8_run/bad_sites", col_names = FALSE)

#now we'd like to generate a whitelist for the populations module of stacks
#this should consist of the snps in th 0.6 version of the vcf excluding the "bad_sites" rows 
#format catalog number \t snp position (from "ID" field)

#create dictionary for converting chrom_pos into stacks catalog ids 
awk 'NR>1870 {print $3}' 0.6.recode.vcf | awk -F :  'BEGIN{OFS="\t"} {print $1,$2-1}' > stacks_loci
awk  'BEGIN{OFS="_";}NR>1870 {print $1,$2}' 0.6.recode.vcf > vcf_loci
paste --delimiters='\t' vcf_loci stacks_loci > catalog_dict

#now remove the bad loci and produce the whitelist
awk 'BEGIN{OFS="\t"} FNR==NR { a[$1]; next } !($1 in a){print $2,$3}' bad_sites catalog_dict > whitelist.txt

7 final populations run

Now that this is all done, run populations using the whitelist

#!/bin/bash

# set max wall-clock time (D-HH:MM:SS)
#SBATCH --time=6-23:59:00

#SBATCH --cpus-per-task=38

# set partition/queue to use
#SBATCH --partition=week-long-cpu

# set name of job
#SBATCH --job-name=popfiltfirst12-8

# set name of output file
#SBATCH --output=popfiltfinalfirst12-8.out

# mail alert at start, end and abortion of execution
#SBATCH --mail-type=ALL

# send mail to this address
#SBATCH --mail-user=ddayan@clarku.edu

source /opt/stacks/2.3/bin/source_me
/opt/stacks/2.3/bin/populations --in-path ./ -M ../info/12-8_good_inds_popmap.txt -t 38  --vcf --merge_sites -e pstI  --whitelist ./filtering/whitelist.txt -H -R 0.95 --ordered-export -O ./final_12-8_filt 

Removed 812327 loci that did not pass sample/population constraints from 962035 loci. Kept 149708 loci, composed of 131008964 sites; 118603 of those sites were filtered, 738822 variant sites remained. 130165888 genomic sites, of which 838483 were covered by multiple loci (0.6%). Mean genotyped sites per locus: 867.70bp (stderr 0.27).

Population summary statistics (more detail in populations.sumstats_summary.tsv): bb: 88.266 samples per locus; pi: 0.057044; all/variant/polymorphic sites: 129901635/738822/317553; private alleles: 68963 cl: 92.487 samples per locus; pi: 0.052766; all/variant/polymorphic sites: 129901635/738822/261896; private alleles: 57993 lb: 77.921 samples per locus; pi: 0.058716; all/variant/polymorphic sites: 129901635/738822/219889; private alleles: 19711 rs: 116.1 samples per locus; pi: 0.061023; all/variant/polymorphic sites: 129901635/738822/560073; private alleles: 281883

So it looks like the final R95 filter removed an additional 131821 SNPs.

Let’s do a quick coverage check and decide on if this is dense enough ## Final coverage check

sampled down to 1000 random sites and ran coverage cehck batch file

gdepths <- read_tsv("log_files/12-8_run/gdepth0.7")
idepth <- read_tsv("log_files/12-8_run/idepth0.7")
ldepthmean <- read_tsv("log_files/12-8_run/ldepth.mean0.7")
gs <- as.data.frame(unlist(gdepths[c(3:370)]))
colnames(gs) <- "depth"

a <- ggplot(data=idepth)+geom_density(aes(x=MEAN_DEPTH))+geom_vline(aes(xintercept = median(MEAN_DEPTH)))+ggtitle("Mean Depth per site per Individual")+xlim(0,50)
b <- ggplot(data=ldepthmean)+geom_density(aes(x=MEAN_DEPTH))+geom_vline(aes(xintercept = median(MEAN_DEPTH)))+ggtitle("Mean Depth per Locus")+scale_x_continuous(breaks = scales::pretty_breaks(n = 10))
c <- ggplot()+geom_density(aes(x=gs$depth))+scale_x_continuous(breaks = scales::pretty_breaks(n = 10))+xlim(0,80)+ggtitle("Depth Per site per individual")

plot_grid(a,b,c)

het miscall rate:

vcf <- read.vcfR(file = "log_files/12-8_run/0.7.sampled.vcf")

## Scanning file to determine attributes.
## File attributes:
##   meta lines: 14
##   header_line: 15
##   variant count: 1000
##   column count: 387
## 
Meta line 14 read in.
## All meta lines processed.
## gt matrix initialized.
## Character matrix gt created.
##   Character matrix gt rows: 1000
##   Character matrix gt cols: 387
##   skip: 0
##   nrows: 1000
##   row_num: 0
## 
Processed variant 1000
Processed variant: 1000
## All variants processed

#now split this into populations to avoid wahlund effect
vcf_popsub <- function(vcf, pop, in_pop) {
  ids <- which(pop %in% in_pop)
  ids <- ids
  vcf <- vcf[, c(1,ids)]

  return(vcf)

}
popid <- unlist(substr(colnames(vcf@gt), 1,2))
vcf_cl <- vcf_popsub(vcf, popid, "cl")
vcf_bb <- vcf_popsub(vcf, popid, "bb")
vcf_rs <- vcf_popsub(vcf, popid, "rs")
vcf_lb<- vcf_popsub(vcf, popid, "lb")
# first get compute expected and observed genotype frequencies
gfreqs_cl <- exp_and_obs_geno_freqs(vcf_cl)
gfreqs_bb <- exp_and_obs_geno_freqs(vcf_bb)
gfreqs_rs <- exp_and_obs_geno_freqs(vcf_rs)
gfreqs_lb <- exp_and_obs_geno_freqs(vcf_lb)

# then plot those.  Set max_plot_loci so that all 2000
# loci will be plotted
geno_freqs_scatter(gfreqs_cl, max_plot_loci = 1000)

geno_freqs_scatter(gfreqs_bb, max_plot_loci = 1000)

geno_freqs_scatter(gfreqs_rs, max_plot_loci = 1000)

geno_freqs_scatter(gfreqs_lb, max_plot_loci = 1000)

#infer het miscall rate
overall_cl <- infer_m(vcf_cl, minBin = 1e15)
overall_cl$m_posteriors

overall_bb <- infer_m(vcf_bb, minBin = 1e15)
overall_bb$m_posteriors

overall_rs <- infer_m(vcf_rs, minBin = 1e15)
overall_rs$m_posteriors

overall_lb <- infer_m(vcf_lb, minBin = 1e15)
overall_lb$m_posteriors

about a 2% miscall rate for heterozygotes